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Abstract 

In this paper we present the study of the math- 
ematical model of a real life joint used in an 
underwater robotic fish. Fluid-structure inter- 
action is utterly simplified and the motion of 
the joint is approximated by Duffing's equa- 
tion. We compare the quality of analytical har- 
monic solutions previously reported, with the 
input-output relation obtained via truncated 
Volterra series expansion. Comparisons show 
a trade-off between accuracy and flexibility of 
the methods. The methods are discussed in 
detail in order to facilitate reproduction of our 
results. The approach presented herein can be 
used to verify results in nonlinear resonance 
applications and in the design of bio-inspired 
compliant robots that exploit passive proper- 
ties of their dynamics. We focus on the poten- 
tial use of this type of joint for energy extrac- 
tion from environmental sources, in this case 
a Karman vortex street shed by an obstacle 
in a flow. Open challenges and questions are 
mentioned throughout the document. 



1 Introduction 

How much of the diverse behavior we observe 
in animals is a direct expression of the dynam- 
ics of the individual's body? In the last two 
decades, many characteristics of animal loco- 
motion on land were successfully linked to the 



mechanical properties of the legs. The springy 
behavior observed in the trajectories of the 
center of mass during running, walking and 
jumping can be explained by the mechanical 
properties of the limbs and its tunning (Dick- 
inson et al. (2000); Farley (1996); Farley et al. 
(1993); Ferris et al. (1998); Kerdok et al. 
(2002); Moritz and Farley (2003); Roberts and 
Azizi (2011)). The hypothesis that running, 
walking and jumping is tuned to the reso- 
nance of the underlying mechanical system is 
strongly supported by experimental evidence 
and by machines constructed based on this 
idea, the passive dynamic walkers (Ahlborn 
and Blake (2002); Alexander (1990); Collins 
et al. (2005); McGeer (1990); Thompson and 
Raiber (1989)). The consequence of such set- 
ting is energy efficient performance and allevi- 
ation of the controller. 

Due to similarities between running and 
swimming (Bejan and Marden (2006); Kok- 
shenev (2010)), it is not surprising (but not 
less exiting) that efficient locomotion in flu- 
ids was reported to relay on the dynamics of 
the body of the animal. Living trouts have 
been observed to exploit the energy in the 
flow they inhabit to reduce their swimming ef- 
forts (Liao et al. (2003)). Later, euthanized 
trouts performed passive self-propulsion when 
placed in the von Karman vortex street shed 
by an obstacle in a flow (Liao (2007)). The 
experimental results are supported by mathe- 
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matical models, analytical and numerical (Al- 
ben (2009); Eldredge and Pisani (2008); Kanso 
and Newton (2009)). These models do not 
fully agree with each other, however this is 
expected due to the mathematical complex- 
ity of the interaction between structures and 
fluids, of which the passive case is the worst 
scenario: unprescribed motion of the interface 
boundary. This alone represents an open chal- 
lenge for the mathematical modeling commu- 
nity. Lest the challenge remains unsolved for 
too long, robotic researchers design their ma- 
chines with less fluid-dynamic rigor, pursuing 
the first passive dynamic swimmer. The final 
objective is to build a swimming machine that 
can perform at least as well as fish, and at com- 
parable power ratings (Harper et al. (1997); 
Lauder et al. (2007)). 

For a robot to extract the energy in the sur- 
roundings, its mechanical properties have to 
be tuned to the environmental energy storage. 
Stated this way, the problem is one of energy 
harvesting, were nonlinear properties are be- 
lieved to be beneficial (Cottone et al. (2009)). 
Therefore, we need to pin down the resonance 
characteristic of the actuators and joints to be 
used in the robot that, as their biological coun- 
terparts, are generally nonlinear. As if diffi- 
culty was lacking, the study of resonance of 
nonlinear systems has suffer a very slow de- 
velopmental process that started early in the 
1960's and has yet not overcome its infancy. 
In a technical report from 1958 by Brilliant 
(1958) we read: 

Sometimes nonlinearity is avoided, 
not because it would have an unde- 
sired effect in practice, but simply be- 
cause its effect cannot be computed. 

Nowadays the situation is not completely dif- 
ferent. However, we have more powerful com- 
puters, new simulation methods and some 
novel uses of classical tools promise to open the 



path ahead Peng et al. (2007); Vakakis et al. 
(2009). 

In this paper we present the study of the 
joint of a robot fish from a mathematical point 
of view. The amplitude of the oscillations of 
the joint in response to periodic forcing is stud- 
ied. In section 2 we briefly introduce the real 
mechanical device and we move to its mathe- 
matical model in section 3. In the same sec- 
tion the model for fluid-structure interaction 
is briefly described. Approximated solutions 
methods are introduced in section 4. Results 
are presented in section 5. Finally, we close 
the paper in section 6 with a discussion on the 
implications of the results and the relevance of 
the approach to the design of robots and the 
test of controllers based on resonances. 

2 A simple compliant joint 

To extract energy from the environment, the 
robotic platform has to be optimally driven by 
interaction forces. In this case, by the interac- 
tion between the rigid body of the robot and 
the surrounding turbulent flow. The feat can 
not be performed if the angle trajectories of 
the joints are fully prescribed by the controller. 
Therefore, the joints of our robot are compli- 
ant, i.e. the motion of the joint is not only 
defined by the controller, but also by exter- 
nal actions. In Figure la we see the details 
of the joints of the robot fish used in Ziegler 
et al. (2011). Each joint behaves rota- 
tional spring. The restoring torque is gener- 
ated when the relative angle between the two 
connected bodies is not zero. The force pro- 
ducing the torque is given by the extension of 
a linear spring fixed to the first body. The 
spring is connected via an inelastic thread to 
an appendage of the second body. When the 
deflection angle is zero, the extension of the 
spring is minimum as well as the force it ex- 
erts. We call tension to this minimum force 
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value and it is referred with the letter F. Mea- 
surements of the torque for F = 0.73 N, are 
given in Figure lb. The values of the param- 
eters used throughout this paper are given in 
Table 1. The two parameters r, d are distances 
that can be seen in the figures. The elastic con- 
stant of the linear spring is K and I denotes 
the moment of inertia of the joint around the 
axis of rotation. The linear specific damping 
coefficient of the joint is denoted with Q. The 
parameters Qo and Q correspond to the am- 
plitude and frequency of the external forcing, 
respectively. 



Name 


Value 




r 


20.24 ± 0.02 mm 




d 


27.68 ±0.02 mm 




K 


81 ± 1 N/m 




I 


(3.1 ±0.1) x 10~ 5 


kg • m 2 


C-i 


(2.2 ±0.1) x 10~ 4 


N • m • s 


Qo-I 


1 x 10~ 4 N-m 






(0,3] Hz 





Table 1: Value of the parameters used here and in 
the model studied in Ziegler et al. (2011) . 



3 Mathematical model 

A geometrical representation of the joint is 
given in Figure 2. The parameters r,d and 
K are fixed at construction time and always 
d > r (Table 1). The tension in the spring 
at its shortest length (F), is the controlled pa- 
rameter and a servomotor can change it dy- 
namically. 

The torque applied to the bodies connected 
to the joint is thus, 



K [ ,/<--- l.Vsin 2 f -e) +F (J) 

-S sm9. 



Where we have defined the parameters e = d—r 
and S = rd. The formula is obtained by calcu- 
lating the deformation of the spring as a func- 
tion of the deflection angle of the appendage 
(9). The deformation is inside parenthesis in 
the numerator of equation (1) and when multi- 
plied by the stiffness K, it gives the force due to 
the deformation of the linear spring. The outer 
factors come from the product of the force and 
the moment arm. Note that the the torque r 
is linear in the controlled input F. 

For 9 <C 1 the third order Taylor expansion 
gives 



t(9, F) 
Where 
k(F) 
a(F) 



K(F)6 + a(F)9 3 + 0(9 5 



j 

e 
e 



2e 2 \F 



(2) 

(3) 
(4) 



And the equation of motion of the deflection 
angle is 



+ C0 + 



9 + C0 + k{F)9 + a(F)9 3 



(5) 



e 2 ±45 sin 2 



Where I denotes the moment of inertia around 
the axis of rotation. The specific damping 
is given by £, and we have defined k = K /i } 
a = a /i. r is the specific net effect of all 
other external torques acting on the joint. The 
approximating equation of motion is the well 
studied Duffing's equation. 

To quantify the error introduced by the ap- 
proximation, we calculated the angle at which 
the difference between the torque produced by 
equation (1) and equation (2) is equal to a 
reference error given by At = rAF, where 
AF = 0.05 N is a reasonable resolution for a 
force sensor working in a 10 N range. These 
angles are plotted in Figure 3 for different val- 
ues of the tension. 
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Figure 1: Description of the joints in WandaX. a) Details of the joints used in Ziegler et al. (2011), 
note the axis of rotation and the appendage, b) Measured torques applied to the joint under controlled 
deflection. Replacing the parameters given in Table 1 into Eq. (1) the fit provides F — 0.73 ± 0.05 N. 




Figure 2: Schematic of the rotational spring used 
to derive Eq. (1). 



3.1 Hardening, linear and softening 
spring 

As can be seen from equation (2), a(F) modu- 
lates the intensity of the cubic non-linear term. 
This term vanishes when 



F = F* 



eK 



(6) 



3S 



+ 1 



rendering a linear spring for the angles where 
the third order approximation is valid. For 
bigger values of F the spring will be softening 
(a < 0) and for smaller values it will be a hard- 
ening spring (a > 0) (Fig. 4). However, the 
full expression of the torque (Eq. (1)) contains 
higher order terms. Hence, the linear behavior 



will be even more evident if the higher order 
terms cancel each other. In Figure 4 curves 
of torque versus angle for several values of F 
are shown. In particular we show the curve for 
which up to seventh order nonlinearities give 
a minimum contribution (found via optimiza- 
tion Fq = 0.9^) together with the curve at 
F = F*. These illustrates the power of the ac- 
tuation chose, since we can control the dynam- 
ical properties of a virtual rotational spring at 
the joint (more details in Ziegler et al. (2011)). 

3.2 Forcing model 

As explained before, the external torques act- 
ing on the joint are due to fluid-structure in- 
teraction. This kind of interaction still poses 
a great challenge for the mathematical model- 
ing community. In the search for simplified 
models of the forces generated in the inter- 
action between flexible bodies and turbulent 
flows, we found the work of Kanso and New- 
ton (2009) and Alben (2009) instructive. Using 
equation (3.8) given by Alben (2009), we can 
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Tension [N] 

Figure 3: Approximation error. Values of the de- 
flection of the joint for which the error reaches 
the reference value given in the text. The 
maximum deflection has a peak close to F3 — 
(5.90 ±0.01) x 10 _1 N meaning that terms of de- 
gree > 3 almost cancel each other. 



-10 1 

6[°) 

Figure 4: Plots of the torque function for values of 
tension F = [0.5, 1, 0.9, 1.5] ■ F*. Note the hard- 
ening behavior for small tension and softening be- 
havior for higher tension. Almost linear curves are 
found for F = F* and F = Fq, as a reference the 
linear curve is shown. 



calculate the pressure difference on the bound- 
ary of a slender body in a vortex street. If 
the width of the body is bigger than the sep- 
aration among vortices and it is placed in the 
middle of the vortex street, the forces on the 
body can be approximated by a sine function 
f(t) = /o sin(fii + <p), where t is time and <j) is 
an initial offset that sets the balance between 
the sine and cosine behavior of the forcing. 
The frequency Q, is proportional to the speed 
of the flow plus the vorticity of the vortices 
(assumed to be equal for each vortex) scaled 
by a factor that depends on the geometrical 
properties of the wake. The amplitude /o is 
proportional to the density of the fluid and the 
square of the vorticity of each vortex. With a 
few additional assumptions, the torque acting 
on the joint can be made proportional to this 
force. Here we adopt this over-simplified forc- 
ing model to avoid diverting the attention of 
the reader from the core ideas of our work. In 
this manner we postpone a detailed study with 
a more elaborated forcing model. 



4 Solution methods 

In this section we present two independent 
methods to estimate the amplitude of oscilla- 
tions of the joint under periodic forcing. The 
first method turns out to be excellent for es- 
timating the amplitude of the first harmonic. 
The second method is based on a Volterra 
series expansion of the Duffing equation. It 
allows to study the response of the joint to 
more general forcing conditions, but the ten- 
sion range for which it is valid is smaller. 

4.1 Harmonic solutions of Dufnng's 
equation 

Under periodic forcing T = Qosin(ilt), equa- 
tion (5) has been extensively studied (see 
Holmes and Rand (1976); Luo and Han (1997) 
and references therein). Following these anal- 
yses, we show here how to maximize the am- 
plitude of the periodic response of the joint by 
tunning the tension parameter. 

The key point of the analysis in Luo and 
Han (1997), is that we search for the ampli- 
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tude of solutions of (5) that are periodic (this 
rules out sub-harmonics, supra-harmonics and 
chaotic motion), Q(t) = Asva(£lt + t/j). Under 
this assumption it can be shown that the am- 
plitude A of such solutions is given by the roots 
of the polynomial, 

9a 2 2 3 (fc-n 2 )a 
16 2 (7) 

A [(ft 2 + c 2 - 2k) n 2 + k 2 } - Ql = 0. 

The roots of this polynomial can be obtained 
analytically using a computer algebra system 
as Maxima (2009) and they establish the re- 
lation between the amplitude of the oscil- 
lations and the parameters of the equation, 
A(k,a,Cl,Qo). In the case at hand, we have 
k(F) and a(F), therefore A(F,n,Q ). For a 
given value of the parameters, only one root 
corresponds to the observed amplitude (there 
are unstable amplitudes). This implies, that is 
not enough to look at the roots, but we must 
also check their stability. This adds some com- 
plexity to the evaluation of the results that 
adds to the limitation of pure harmonic inputs. 

4.2 Volterra series expansion 

When the amplitude of the forcing is suffi- 
ciently small, the behavior of Duffing's oscil- 
lator in the neighborhood of the origin can 
be described by polynomial Volterra function- 
als(Theorem 3.1 of Rugh (1981)). This means 
that the angle of the oscillations can be de- 
scribed by an expansion of the form 

n 

9(i)»E»0- ( 8 ) 

i=l 

where n is the order of the expansion and each 
term is the multi-dimensional integral 

/+oo 
hi(t - (71, . . . ,t- Oi) 
-oo 

T{a 1 )...T(a l )da 1 ...da i . (9) 



where T(t) is the forcing signal and we wish to 
determine the kernel functions {/i«}™ =1 - How- 
ever, for the objective at hand, it is more use- 
ful to calculate the Fourier transform Y(u>) 
of these functionals. We have not found this 
process described in the literature, hence we 
describe it succinctly. We used the inverse 
Fourier transform definition on the input sig- 
nal and substitute it into the functional defi- 
nition Eq. (9). By inverting the order of in- 
tegration and splitting variables, we make the 
kernel transform H{uj) appear. Then, we com- 
pare this expression to the formal definition of 
the inverse Fourier transform of Y(cj) and iso- 
late the desired result. This process can be 
applied for the first three orders, but it gets 
cumbersome for higher ones. Therefore, be- 
fore showing the results, we will briefly intro- 
duce some notation used here to simplify the 
presentation. 

An ordered set of arguments (si, S2, ■ ■ ■ , s n ) 
will be denoted si :n . In general we have, 

(s k ,s k+1 ,. . . ,s n ) = s k:n k<n. 

For example, the fifth order kernel 
Hs(sx, S2, S3, S4, S5) will be written 

#5(S1:5) = 

- 3a #i (Ss 1:5 ) ffi(*i)ffi(s 2 )ff 3 (s3:5), 

where Ssi : 5 = Yli=i s i- For the integration 
variable in multiple integrals we will write 
dsi:n = dsids2 • • • ds n . For multiple summa- 
tions over n indexes we will write Xyjt . = 
Sfc x • • • Ylk n - Having defined the notation, we 
continue with our presentation. 

For zero initial conditions, (6, 6) = 0, and 
working in the frequency domain, each kernel 
can be obtained recursively based on lower or- 
der kernels (details of the calculation are given 
in Peyton- Jones and Billings (1989)). The first 
order kernel derived from the model in Eq. (5), 
is equivalent to the frequency response func- 
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tion (also known as transfer function) of a lin- 
ear second order system, 

*<•> = (10> 

The following equations present the Volterra 
kernels up to seventh order using recurrent re- 
lations, 

#3(S1:3) = -aHi (Ssi:3) ■ 

H5(sx:5) = SaHi (Ssi :5 ) • 
Hi(s\:7) = - 3aHi (Ssi :7 ) • 

H 1 {s l )H 1 {s 2 )H 5 { S3:7 ) + 
H 1 (s 1 )H 3 (s 2 . A )H 3 {s 5 .. 7 ) • 

Due to the symmetry of the eq. (5), all even 
order kernels are zero. Note that these kernels 
are valid for any second order cubic oscillator. 
The recursive formulas were obtained using the 
probing method of Peyton-Jones and Billings 
(1989), which proceeds as follows. We start 
from a polynomial nonlinear ordinary differen- 
tial equation (as Eq. (5)) and assume the re- 
sponse can be represented by Volterra series, 
i.e. Eq. (8). We substitute this into the sys- 
tem's differential equation, and equate similar 
terms. This last step produces equations relat- 
ing different order kernels and the expressions 
are always recursive (we obtain Eqs. (11) in 
our case). The factorial in the recursive rela- 
tions appear from the multinomial expansion 
of nonlinear polynomial terms. 

We proceed to calculate the Fourier trans- 
form of the output given by Eq. (9) subject to 
harmonic forcing. We will use Eqs. (11) in or- 
der to obtain an explicit relation between the 
parameters of our model k, a, (, Qq, Q, F and 
the amplitude of the first harmonic of the out- 
put. 

We recall that the Fourier transform of a 
general bandwith limited periodic function is 



a scaled Dirac Comb (Shah function, Impulse 
train, Dirac train, etc.) (see Schetzen (1980)), 

+L 

= ^ XiS(u-m) (12) 

l=-L 

where L is a positive integer, is the complex 
amplitude of the i-th harmonic and <5( • ) is the 
impulse function (Dirac delta function). The 
Fourier transform of each of the output terms 
in (8) is 

/oo 
Hi (Awi) • 
-oo 

y(wi) ■ ■ ■ 7(^-1)7 (u) - Ewi-.i-i) du>i :i _i = 

+L 

(2tt) 1 - < Yl B h:h 5 (u - nU 1:i ) , 
h-.i=—L 

(13) 

where we have used Eq. (12) and Aa;j = 
(wi, . . .,Wi-i,uj - Ewui-i) and, 

B h:l .=H i (M 1:i )X h ---X l .. (14) 

The integral in (13), has an interesting geo- 
metric interpretation in terms of convolutions 
in hyperplanes, we refer the interested reader 
to Lang and Billings (1996). Therein, the out- 
put frequency range of nonlinear systems that 
are representable by Volterra series is analyt- 
ically calculated. Additionally, in Eq.(16) of 
that paper the frequency spectrum of the out- 
put signal is represented as the superposition 
of contributions from the nonlinearities. In 
the case studied herein, the input has only 
one frequency, therefore L = 1, Xq = 0, 
X±i = ±jQoir with j the imaginary unit. Re- 
placing these values in all the equations and 
using the relations in Eq. (11) we obtain the 
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Figure 5: Amplitude of periodic response in the plane (F, f2). (a) Amplitude of the oscillations according 
to Eq. (7). The line of maxima is marked with dashes. The vertical dotted lines indicate selected values 
of Q. (I)) Relative differences between the predicted amplitudes by the method described in section 4.1 
and the Volterra expansion described in section 4.2. 



desired result, 
m = ^{3a 3 Q 7 [ 

+ 6H 1 (3n)Hf{-n)Hf(n) 
+ 6Hf(3n)Hf(-n)Hf(n) 
+ 3.ffi(-3n)iyf(-n).Hf(«) 

- 15.ffi(3fi).H?(-«).Hf(fi) 

+45Hf(-n)Hf(n)+18Hf(-n)Hf(fl) ] 

- 12a 2 Q 5 [Hx(3n)Hf(-n)Hf{n) 

+ 6H^(-n)Hf(n) 

+ 3Hf(-n)Hf(Q)] 

+ 48aQ^i(-^)if 1 3 (0) - 64Q #i(^)}- 

(15) 



Where H\{x) is given by Eq. (10) with s = jx. 
The amplitude of the oscillation is obtanied by 
taking the double of the modulus of the com- 
plex number Y(U)/tt. 

5 Results 

Figure 5a shows the amplitude of periodic os- 
cillations in the (F, 17) plane, according to 
Eq. (7). The line of maxima is shown with 
dashes. This line describes the value of the 
tension that produces maximum amplitude for 
a given forcing frequency. In Fig. 5b we plot 
the logarithm of the relative difference between 
the amplitudes given by (7) and (15). The 
two approximation differ for regions of low fre- 
quency and low tension where the system is 
most nonlinear (Luo and Han (1997)). To com- 
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pare with the observed amplitude (obtained by 
simulating (5) without any approximation), we 
extracted the curves of amplitude against ten- 
sion for ^/2tt = 0.5, 1.5 Hz (vertical dotted lines 
in Fig. 5). These curves are shown in Figure 6 
together with the amplitude calculated using 
Eq. (15) corresponding to the Volterra ker- 
nel expansion. For the 1.5 Hz frequency, both 
models predict the simulated amplitude accu- 
rately. For the 0.5 Hz frequency, the ampli- 
tude predicted by Eq. (7) drifts away from the 
simulated value for lower tensions. The ampli- 
tude calculated using the Volterra expansion 
diverges for low tensions at this frequency. In 




2 4 6 8 10 



F[N] x lO- 2 

Figure 6: Amplitude of periodic response of the 
joint for forcing with frequencies O = [0.5, 1.5 ]Hz. 
Amplitudes according to Eq. (7) in solid line, 
amplitudes obtained from the Volterra expansion 
Eq. (15) in dashes and amplitudes from simulations 
without approximation in circles. 

the same figure we show the amplitudes mea- 
sured in simulations using the exact expression 
for the torque (1) (Octave (Eaton (2002)) func- 
tion ode45, absolute and relative tolerances, 
10 -6 ). The amplitude of the oscillations is ob- 
tained using the absolute value of the analytic 
signal of the angle (Octave function hilbert). 
The agreement between the approximated re- 
sults and the simulated ones is noteworthy. 
These results show that the approximation 



from Luo and Han (1997) can be used for a 
tension controller designed for a joint of this 
kind, to obtain periodic responses with the 
same frequency as the forcing. The Volterra 
approximation fails when the system becomes 
strongly nonlinear, however these expansion 
can be used for more general responses or when 
the inputs have a more complicated frequency 
spectrum (as in a realistic scenario). 

6 Discussions and conclusion 

Though the quote from Brilliant (1958) re- 
mains valid, we used the knowledge about 
Dumng's equation to understand analytically 
a compliant joint. This gives a corner stone 
to study more complicated setups (e.g. chain 
of joints, as in the robotic fish) and define the 
robust engineering design of robots with the 
desired resonance properties. Additionally, the 
results showed here provide a reference solution 
to the problem of finding the right pretension 
for a given external forcing, related to the idea 
of adaptive controllers. For a forcing with a 
sufficiently slow varying frequency, the tension 
could be adjusted to maintain the amount of 
energy transferred into the system to its maxi- 
mum possible value, i.e. keep the system close 
to the line of maxima in Fig. 5a. With the 
results herein the adjustment could be done 
by direct calculations. However, other meth- 
ods like the frequency oscillators(AFO) (Buchli 
et al. (2006)) were used to achieve a similar 
objective. Applications of the latter has been 
reported in Buchli and Ijspeert (2008), how- 
ever, due to lack of theoretical results on the 
resonance frequency of the nonlinear platform 
studied therein, the theoretical solution to the 
problem was unknown and a grounded evalua- 
tion of the performance was not possible. Ap- 
proximations as the one presented here could 
be used to verify the results generated by the 
adaptive oscillators. 
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In our results we exclude any information 
about oscillations with frequencies different 
from that of the forcing since that is the case 
where Eq. (7) is valid. The higher order har- 
monics in the system's output could be deter- 
mined and compared with the simulated re- 
sults. Traditional harmonic balance method 
can not be used to find these results. Ad- 
ditionally, a more complete panorama of the 
frequency response of the system could be ob- 
tained with a nonlinear analysis method as the 
NOFRF (proposed in Lang and Billings (2005) 
and showcased in Peng et al. (2007)), which is 
based on Volterra series expansions as the one 
presented in this work. NOFRF give the re- 
sponse of the system to a given set of inputs, 
therefore the natural ensemble of forcing sig- 
nals must be known to acquire useful informa- 
tion. This ensemble has to be compiled from 
(or modeled based on) fluid-structure interac- 
tions data. Comparing the present results with 
the information provided by NOFRF will be 
the objective of our next work. Proved that 
such analysis is helpful for the design of com- 
pliant robots, we will extend the analysis to 
biped and quadruped robots. In that situa- 
tion the forcing comes from the ground reac- 
tion forces, that are generally not smooth (im- 
pacts) and contain a rich frequency spectrum. 
We consider those scenarios as highly complex 
problems to be addressed with a mature tool- 
box of methods. 

Concluding, we have provided a method that 
allows the construction of a controller that 
would maximize the harvesting of environmen- 
tal energy, under the circumstances defined 
herein. The method will be extended to cover 
more realistic situations than the simplified 
periodic forcing, including forcing with wider 
frequency spectrum. Additionally, the same 
methodology presented here can be used to 
verify results obtained using heuristic methods 
as AFO or neural networks. Our results can be 
easily validated and we invite our colleagues to 



do it. 
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